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While evaporating solvent is a widely used technique to assemble nano-sized objects into desired 
superstructures, there has been limited work on how the assembled structures are affected by the 
physical aspects of the process. We present large scale molecular dynamics simulations of the 
evaporation-induced assembly of nanoparticles suspended in a liquid that evaporates in a controlled 
fashion. The quality of the nanoparticle crystal formed just below the liquid/vapor interface is found 
to be better at relatively slower evaporation rates, as less defects and grain boundaries appear. This 
trend is understood as the result of the competition between the accumulation and diffusion times of 
nanoparticles at the liquid/vapor interface. When the former is smaller, nanoparticles are deposited 
so fast at the interface that they do not have sufficient time to arrange through diffusion, which 
leads to the prevalence of defects and grain boundaries. Our results have important implications in 
understanding assembly of nanoparticles and colloids in non-equilibrium liquid environments. 



I. INTRODUCTION 

Large, defect-free crystalline arrays of nanoparticles 
(NPs) are critical for many technologically important 
ultra thin film materials including filters, sensors, opti- 
cal devices, and magnetic storage media. 1 2 One common 
method for dispersing nanoparticles is to suspend them 
in solution, spread the suspension on a surface, and then 
evaporate the solvent SEMI This technique has been used 
to assemble NPs into desired structures such as nanoclus- 
ters, rings, wires, stripes, films, and super latt ices PMS1 
The factors that influence the assembled structure in- 
clude the evaporation kinetics of the solvent, the flow 
properties and concentration of NPs in the solution, the 
NP-NP interactions, the interactions between NPs and 
the liquid/ vapor interface, and the wetting/dewetting be- 
havior of the suspension on the solid surface. 15 One ubiq- 
uitous example is the " coffee-ring" stain left behind after 
a particle-laden drop has evaporated. A simple and beau- 
tiful physical picture of this class of phenomena was first 
clarified by Deegan et alP^ When the drop dries, a flow 
of liquid from the interior of the drop to its edge emerges 
to replenish the liquid that has evaporated from the re- 
gion around the pinned contact line. This flow convects 
particles to the drop's edge, where they form ring-like de- 
positions. The model was tested recently by Chen et al. 
on the molecular scale with molecular dynamics (MD) 
simulations.^ 

Recently, Marin et al. demonstrated an order-to- 
disorder tr ansiti on of the particle organization in coffee- 
ring stainsPH^ In the early stage of evaporation, the 
deposition speed of particles at the contact line is low 
and they have time to arrange into an ordered super- 
lattice by Brownian diffusion. However, at the end of 
evaporation, because of the increase of the flux velocity 
toward the contact line, particles are deposited at a high 
speed and jammed into a quenched disordered phase as 
the result of "rush-hour" traffic. The experiment clearly 
showed the importance of deposition/accumulation and 
diffusive time scales on controlling the particle ordering 
during solvent evaporation. 



Bigioni et al. showed earlier that particle accumulation 
and diffusion times determine the assembly behavior and 
growth kinetics of monolayers of gold NPs w ith diame- 
ters ~ 6 nm at a liquid/vapor interface! 6 * 10 * 23 * They found 
that in order to induce NP assembly at the interface, the 
evaporation rate of the solvent has to be rapid and an 
attractive interaction between the NPs and interface is 
required. The latter was induced by the excess of do- 
decanethiol in a toluene solution where gold NPs were 
suspended. After the assembly is initiated, the growth 
kinetics of each crystalline domain can be exponential, 
sub-exponential, or linear in time, depending on the ra- 
tio of a diffusive length scale to the domain size. By 
tuning the concentration of dodecanethiol, Bigioni et al. 
were able to make long-range ordered arrays of gold NPs 
over macroscopic areas 

The work of of Bigioni et al. showed that evaporation 
generally needs to be fast so that the accumulation of NPs 
at the interface passes certain critical density, leading to 
the nucleation and growth of 2-dimensional (2D) crystals. 
However, if evaporation is too fast, the crystalline qual- 
ity may diminish because NPs do not have sufficient time 
to arrange via diffusion, but instead quenched disordered 
structures will emerge, similar to the "rush-hour" effect 
observed in the coffee-ring stains. This trend was actu- 
ally observed earlier in an experiment of Im and Park 
using larger polystyrene colloids with diameters ~ 230 
nm j24| w j iere the evaporation rate was tuned by varying 
the temperature. They found that there exists an opti- 
mum evaporation temperature that yields particle arrays 
of the highest quality at the water /air interface. When 
the evaporation temperature is higher than the optimum 
one, the accumulation rate of colloids at the interface is 
too large, leading to defects in the array. It is thus in- 
teresting to determine if the effect of evaporation rate 
on particle assembly persists to smaller NPs in the fast- 
evaporation regime, which nicely fits the accessible length 
and time scales of MD simulations. This was the moti- 
vation of our work reported here. To anticipate the final 
results, we will show that there exits an optimum rate, 
determined by the diffusive times of NPs along the in- 
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terface, at which the 2D crystal of NPs formed near the 
interface has the highest quality even without the "an- 
nealing" step. 

The remainder of this paper is organized as follows. 
The simulation methods are briefly introduced in Sec. II. 
Then Sec. Ill is devoted to results and discussion. Con- 
clusions are included in Sec. IV. 

II. SIMULATION METHODS 

We modeled the liquid solvent and its vapor with 
atoms interacting through a standard Lennard- Jones 

— where 



O/r) - (cr/r) 



(LJ) potential, U^(r) = 4e 

r is the distance between two atoms, e the unit of en- 
ergy, and a the diameter of atoms. The potential was 
truncated and shifted to at r c = 3.0<j. The mass of 
each L J atom is m and a time unit r can be defined as 
r = ^Jma 2 Je. The NPs of diameter d = 20a were as- 
sumed to consist of a uniform distribution of atoms inter- 
acting with a LJ potential. For spherical particles, their 
mutual interaction can be determined analytically by in- 
tegrating over all the interacting LJ atom pairs between 
the two particles) 2 ^ The resulting potential depends on 

the inter-NP distance, and a Hamaker constant A nn 
that characterizes the interaction strength. We used a 
standard value A nn = 39.48e and a cut-off 20.427a. This 
led to a purely repulsive NP-NP interaction, which physi- 
cally corresponds to adding a short surfactant coating on 
NPs to avoid flocculation! 26 * 27 ! The interaction between 
LJ atoms and NPs was determined similarly by integrat- 
ing the interaction between a single LJ atom and those 
within a NP. We set the corresponding Hamaker constant 
A ns = lOOe, which with a cut-off of 14a resulted in NPs 
that were fully solvated and dispersed in the LJ solvent.^ 

The simulation cell was a rectangular box of dimen- 
sions L x x L y x L z . The liquid/vapor interface was par- 
allel to the x-y plane, in which periodic boundary con- 
ditions were employed. In the z direction, the LJ atoms 
and NPs were confined by two flat walls at z = and 
z = L Z: respectively. The depth of the liquid was at 
least 150a and the vapor was at least 70a thick before 
evaporation. The resulting liquid and vapor densities are 
0.64ra/a 3 and 0.056ra/a 3 , respectively. Note that LJ 
monomers have a vapor density higher than most real 
fluids (Fig.JTJa)). This implies that the vaporization rate 
of LJ liquids is quite large and they can evaporate very 
fastP^ To provide a rough mapping of the evaporation 
rates (number fluxes of solvent atoms) in our simula- 
tions to experimental values, we take a ~ 0.3 nm and 
r ~ 2 ps, the slowest evaporation rate we can simulate 
with MD is at the order of 10 _4 r _1 a -2 , corresponding 
to ~ 5 x 10 26 sec -1 m -2 . This is about an order of mag- 
nitude higher than the maximum rate at which water 
evaporates, which is not surprising since water has an 
unusually high boiling temperature and surface tension 
compared with the LJ liquids studied here. 

All simulations were performed using 

LAMMPSPMI 

The equations of motion were integrated using a velocity- 



Verlet algorithm with a time step St = 0.005r. During 
the equilibration, the temperature T was held at 1.0e/&B 
by a Langevin thermostat weakly coupled to all LJ atoms 
with a damping constant T — O.lr -1 . Once the liq- 
uid/vapor interface was equilibrated, the Langevin ther- 
mostat was removed except for those liquid atoms within 
15a of the lower wall at z = 0. We confirmed that the 
evaporation occurring at the interface was not affected 
by the thermostat.^ 

III. RESULTS AND DISCUSSION 

A. Effect of Evaporation Rate 

Most of our simulations started with 668 NPs dispersed 
in a liquid that contained about 7.1 million LJ atomsj^ 
as shown in Fig. [TJa). The volume fraction of NPs was 

0. 205, which was based on the diameter of bare NPs. 
If we took into account the effect of a layer of solvent 
of thickness ~ la that was effectively attached to NPsJ^ 
the effective volume fraction of NPs would be even higher, 
roughly by a factor of 1.16. We have studied systems 
of which the liquid/vapor interface has various aspect 
ratios. Since the NPs were spherical and with a short- 
range isotropic interaction, their ordering at the inter- 
face was expected to be hexagonal. Thus to reduce the 
finite size effects, we focused on a system with a cross- 
section that was commensurate with hexagonal packing, 

1. e., L x = VSL y /2 = 271a (Figs, [j] and 2l. L x and L v 
were chosen so that the lattice spacing in a perfectly or- 
dered hexagonal array would be the same as that found 
after some of the solvent was evaporated in a system that 
had a square cross-section with L x = L y = 304a. How- 
ever, this latter system has a larger interfacial area than 
the previous one. To reduce the potential influence of 
areal size, we studied another system of which L y was 
increased to 352a while L x was kept at 271a, so that not 
only the interfacial area was close to that of the square 
interface, but also the aspect ratio of the interface fa- 
vored a hexagonal packing with a lattice spacing found 
in the previous two systems. Our results indicated that 
the size and aspect ratio of the interface had negligible 
effects on the particle packing. Below we mainly discuss 
the results for the first system and include the results of 
other two systems in Supplementary Material.^ 

The side view of the NP distribution in the equilibrium 
solution is shown in Fig. [TJe), where the dashed line in- 
dicates the location of the liquid/ vapor interface. NPs 
are almost uniformly distributed in the solvent. Note 
that during all our simulations, including the evapora- 
tion runs described below, all NPs were fully solvated in 
the liquid solvent. NPs within d = 20a of the interface 
were grouped together and referred as the top layer. At 
equilibrium their areal coverage is about 45%, which is 
below the critical coverage 70% at the 2D hexatic phase 
transition.^ As shown in Fig. [lji), the layer is clearly 
disordered. 

To quantitatively characterize the in-plane packing ge- 
ometry, we have defined an order parameter iJjq for each 
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FIG. 1: Snapshots of a NP solution with L x = \^3L y /2 = 271 a: (a) the equilibrated solution before the evaporation of 
solvent (containing ~ 7.1 million LJ atoms); (b) after 52% (~ 3.7 million atoms) of the solvent had evaporated at a fixed 
rate 2.35 x 10~ 4 r~ 1 a~ 2 for 1.85 x 10 5 r; (c) after 52% (~ 3.7 million atoms) of the solvent had evaporated into a vacuum for 
2.9 x 10 4 t; (d) the same system as in (c), but the system had re-equilibrated for 3.4 x 10 4 r after the evaporation was stopped. 
For clarity in (a)-(d) only about a quarter of the system in the x-y plane is shown, (e)-(h): Side views of the NP distribution 
for systems in (a)-(d), respectively; the dashed lines indicate the location of the liquid/vapor interface, (i)-(l): Top views of 
NPs in the top layer for systems in (a)-(d), respectively; colors are based on an order parameter (ipe) for each NP as defined 
in Sec. III(A). 



NP. A Voronoi construction was performed and those 
NPs with 6 neighbors were identified since in a close 
2D packing each particle was expected to have 6 nearest 
neighbors forming a hexagon. The angle Oi correspond- 
ing to the i-th side of the hexagon was found and iJjq was 
computed as ipQ = Y^=i 1^ — tt/6| . We found that ipQ 
was typically in the range of to 4tt/3. So we normalized 
by 47r/3 and the normalized was thus in the range 
of to 1. ^6 = indicates a perfect hexagonal packing, 
while ^6^1 indicates situations far from a hexagonal 
lattice. The Voronoi construction also showed that all 
other NPs had either 5 or 7 neighbors and they were 
assigned i/jq = 1 to indicate that their packing configura- 
tion was not hexagonal. In Fig. ITVi)-(l), NPs are colored 
according to their values of i/jq (NPs with i/jq < 0.2 are 
colored red and with ip$ = 1 are colored blue). As shown 
in Fig. [TJi), in an equilibrated solution NPs near the in- 
terface are randomly distributed without any ordering, 
which is also clear from the Voronoi construction shown 
in Fig.^a). 



Evaporation of the solvent was modeled by introduc- 
ing a deletion zone from (L z — 20<j) to L z . All or a 
certain number of vapor atoms in the deletion zone were 
removed every 0.5r. The former mimicked evaporation 
into a vacuum which proceeded at the maximum rate, 
while the latter simulated evaporation at a fixed rate 
2e-> namely the rate of removal of solvent atoms out of 
the simulation box. Results in Figs. []Jb), (f), and (j) 
are for a system of which the solvent had evaporated at 
je — 2.35 x 10 _4 r _1 <j -2 since t = 0. During evapora- 
tion, the liquid film where NPs were dispersed shrunk, 
and NPs started to accumulate below the liquid/vapor 
interface that impinged from above. In this case the va- 
por density was only slightly reduced compared with its 
equilibrium value (Fig. [TJb)) throughout the whole sim- 
ulation and the interface always remained close to equi- 
librium. As a result, the rate je was roughly the same 
as the actual outward flux of solvent atoms near the liq- 
uid/vapor interface or in the vapor. The side view of the 
NP distribution in Fig.[]Jf) shows that NPs formed layers 
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where their concentration had peaks. At t = 1.85 x 10 5 r, 
52% of the solvent (~ 3.7 million atoms) had evaporated, 
and the volume fraction of NPs increased to 0.35. At 
least 3 layers of NPs were formed as shown in Fig. [ljf ) , 
with 171 NPs in the top layer. Figure [TJj) shows the 
in-layer structure of these NPs, from which a hexagonal 
packing is clearly visible. This is also evident from the 
corresponding Voronoi construction (Fig. [2jb)). 

Note that the number of layers induced by evaporation 
depends on the volume fraction of NPs and the interfacial 
area. In many experiments, NP solutions are quite dilute 
and only one layer of NPs forms, but more layers were 
also experimentally observed, 6 as in our simulations. The 
structures in the NP layers below the top one were essen- 
tially similar but with a diminishing crystalline quality 
as the distance from the interface increased. Two con- 
secutive layers were roughly in-registry as in fee or hep 
crystals, but we had too few layers to distinguish the two 
crystal structures. For clarity, we focus on the in-plane 
structure in the top NP layer below. 

The results in Figs, [ljc), (g), and (k) are for a sys- 
tem with the same starting state, but with the solvent 
evaporating into a vacuum since t = 0. Here the va- 
por was quickly depleted (Fig.[TJc)) and the evaporation 
at the liquid/vapor interface proceeded at the maximum 
rate, which decreased with time and eventually reached a 
plateau value around 1.2 x 10 -3 t -1 <t -2 (see Fig{3|. The 
actual flux of solvent atoms at the interface was lower 
than je for the first lOOOr during which the vapor de- 
pletion occurred, and then the two were almost identi- 
cal. It only took 2.9 x 10 4 r to evaporate roughly the 
same amount (52%) of solvent as in the previous case 
with a fixed rate. The layering of NPs was obvious from 
the side view (Fig. [TJg)). However, each layer was com- 
posed of several 2D grains separated by grain boundaries, 
where the packing deviated significantly from a hexago- 
nal lattice as indicated by colors based on the values of ip$ 
(Fig. Qk)) and the corresponding Voronoi construction 
(Fig. |2[c)) for NPs in the top layer. 

Comparing Figs. [TJj) and (k) (or Figs. |5|b) and (c)) 



shows that the slower evaporation led to a NP packing 
with a better crystalline order. The same trend was also 
found for the other two systems evaporating at various 
rates (see Figs. SI (a) and (b), and Figs. S2(a) and (b) 
in Supplementary Material). 33 This effect of evaporation 
rate ]e on the assembly quality can be understood on 
the basis of two time scales: one describes how fast NPs 
accumulate at the interface, and another is determined 
by the diffusivity of NPs along the interface. In the 
steady state, the liquid/vapor interface recedes almost 
uniformly, with a receding velocity v that is related to 
3e through 2e — pv, where p = 0.64a -3 is the liquid 
density. Using a model discussed by Bigioni et. a/.p^a 
time scale t a = (cjyvd 2 ) -1 can be defined to describe the 
speed of NP accumulation at the interface, where <p is the 
NP concentration in the solution. After NPs reach the 
top layer, they primarily diffuse along the interface. We 
denote the diffusive time scale as td, which can be esti- 
mated as td = d 2 /AD with D as the diffusion coefficient. 
For d = 20a and D^3x 10~ 3 cr 2 /V, 28 the diffusive time 
is t d ~ 3.3 x 10 4 r. 

NPs accumulate faster than they diffuse in the limit 
ta <^ td, where the growth of 2D crystals is induced pri- 
marily by NP addition. As it is unlikely that crystalline 
grains in different spatial regions have orientations that 
match, this leads to the occurrence of grain boundaries 
when grains meet. However, in the opposite limit t a > 
NPs are added to the top layer slowly and small crystal 
regions can act as seeds which then grow larger through 
NP diffusion. In this case defects and grain boundaries 
are less likely to appear. The NP layer also has more 
time to relax into a uniform hexagonal lattice even if im- 
perfections in the packing occur. For the system shown 
in Figs. []Jk) and[2|c), the solvent evaporated fast and 
the interface velocity v was about 2 x 10 _3 <j/r, which 
leads to t a ~ 1.3 ~ 2.3 x 10 4 r. The value of t a has a 
range since the NP concentration <fi increased as the sol- 
vent evaporated. For Figs. [TJj) and [5Jb) where ]e was 
fixed at 2.35 x lO^r" 1 ^ 2 , v was about 3.7 x 10- 4 a/r, 
which gives t a ~ 7 ~ 12 x 10 4 r. As expected, t a < td 
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in the fast evaporation (into a vacuum), while t a > in 
the slow evaporation (at a small fixed rate). 

Since t a and ]e are related through t a = (<j)d 2 j e / p)~ x , 
by equating t a to td one can estimate a critical rate 
jc — p/{(j>d 2 td)- When j E < jc, an almost defect- 
free hexagonal NP array is expected; while for j E > jc 
grain boundaries and point defects will appear in the 
2D crystal induced by evaporation. For our simulations, 
jc ^9x1CT 4 t -1 <7 -2 . When the solvent evaporated into 
a vacuum, the evaporation rate was always larger than 
jc, leading to NP layers with defects and grain bound- 
aries. However, when j E was reduced and fixed at a 
value below jc, the assembly quality in NP layers was 
greatly improved. Furthermore, from experiments it is 
already known that if the evaporation is too slow, then 
the density of NPs at the interface never reaches the crit- 
ical nucleation density since NPs have abundant time to 
diffuse back into the bulk of the solution, and assembly 
is not initiated in the first place!^^ Combined with this 
observation, our simulations indicate that there exits an 
optimum evaporation rate j opt at which the NP arrays 
formed near the liquid/ vapor interface are of the highest 
quality. Though it is very difficult to precisely calcu- 
late jopt from simulations, our results indicate that j opt 
should be about an oder of magnitude smaller than jc 
since only a few defects and grain boundaries were ob- 
served at j E — jc/5 (Fig. 2^b) and Fig. Sl(a) in Supple- 
mentary Material), 33 while more were seen at jE — jc 
(Fig. S2(a) in Supplementary Material) P2 Therefore, a 
reasonable estimate is j opt ~ jc/10 = p/ (10 (j)d 2 td). 

B. Effect of Annealing 

Since the volume fraction of NPs was only 0.35 in 
Fig. [TJc) and still below the critical fraction of hard- 
sphere crystallization (~ 0.545)p^we expected that the 
layering structure in Fig. [TJg) would be destroyed when 
the evaporation was stopped and the system was allowed 
to relax. However, this was not the case on MD time 
scales, and the layering became even more dramatic af- 
ter the evaporation ceased at t = 2.9 x 10 4 r and the 
system had re-requilibrated, as shown in Fig. [ljd) and 
(h) at t = 6.3 x 10 4 t. The robustness of layering was 
partly due to the fact that some solvent evaporated to 
fill the vapor-depleted region above the liquid and to re- 
establish the liquid/ vapor equilibrium. This increased 
the volume fraction of NPs in the solution to 0.38. Fur- 
thermore, melting requires the diffusion of NPs in the di- 
rection perpendicular to the layers, which was hindered 
by the ordered structure in each layer and the commen- 
surability between layers. Our results thus indicate that 
the system was dominated by non-equilibrium dynamics 
with long relaxation and equilibration times, much larger 
than those achievable with MD simulations. 

The organization of NPs in the top layer was greatly 
improved after relaxation. The grain boundaries in 
this layer completely disappeared after approximately 
3.4 x 10 4 r and NPs formed an almost perfect hexag- 
onal lattice, as shown in Fig. [TJl) and Fig. j2^d) (also 



see Figs. Sl(c) and S2(c) in Supplementary Material) 
This finding is consistent with an experimental observa- 
tion that alternative evaporation and relaxation produce 
2D crystals of NPs with less defects. 36 It reveals the im- 
portance of relaxation that allows imperfections in the 
NP packing to heal via structural re- arrangement, which 
only gradually proceeds through the diffusion of NPs. 
The improvement of the crystalline quality was also ob- 
served for the other layers below the top one after anneal- 
ing, though over the time scale of simulations they were 
still not completely ordered. For NPs with diameters at 
the nanometer scale and without mutual attractions as 
studied here, their diffusion is relatively fast and thus 
the time scale of structural re-arrangement is accessible 
with MD. However, it can be imagined that for larger 
NPs or colloidal particles in the micrometer range and 
with attractive inter-particle interactions, diffusion be- 
comes much slower and it will take much longer time to 
remove defects and grain boundaries. This is one of the 
reasons that they are frequently observed in assembly 
experiments 

C. Effect on Evaporation Rate of Assembled 
Nanoparticle Layers 

When NPs crystallize below the liquid/vapor interface, 
they slow down the evaporation, since liquid atoms need 
to circumvent NPs to reach the interface and then evap- 
orate. This blockage effect on evaporation is shown in 
Fig. [3^ a), where the rate jE is plotted against time for a 
neat solvent and the NP solution that both evaporated 
into a vacuum. At earlier times, NPs were dilute near 
the interface and liquid atoms essentially did not feel 
their presence. Therefore jE was roughly the same for 
the two systems. At later times, since NPs started to 
form a dense layer below the interface, jE for the NP 
solution became smaller than that for the pure solvent. 
The factor of reduction was 25% when jE reached the 
plateau values in the two systems at t > 800r. How- 
ever, the reduction in jE is not as dramatic as what a 
naive estimate of the areal reduction of the liquid/vapor 
interface would indicate. The reason is that NPs still sit 
below the interface in our simulations and their top sides 
only barely touched the liquid/vapor transition zone. As 
illustrated in Fig. |3jb) , only liquid atoms attached to the 
surface around the north pole of NPs were affected in 
terms of evaporation. At several a away from this region 
(outside the box enclosing the north pole of the immersed 
NP in Fig. |3jb)), the distance from the interface to NP 
surface quickly increases beyond a few diameters of LJ 
atoms, and the evaporation of the liquid there was not 
influenced by the presence of NPs. However, if the NP 
layer straddled the liquid/vapor interface (see the strad- 
dling NP in Fig. [3jb)), we would expect jE to decrease 
by an order of magnitude compared with the case of a 
neat solvent. In this case, the interfacial area is indeed 
reduced by NPs. The evaporation cannot occur in the 
area covered by NPs and will slow down significantly. 
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FIG. 3: (a) Evaporation rate je vs. time for a NP solution 
(solid red line) and a neat solvent (dotted green line), showing 
the blockage effect of NPs on evaporation, (b) A schematic 
illustration of how the evaporation is affected by nanoparticles 
immersed just below or straddling the liquid/vapor interface: 
only the evaporation of the solvent in the area enclosed by 
the black boxes is blocked. 

IV. CONCLUSIONS 

To summarize, we have studied the assembly of NPs in 
a solution induced by solvent evaporation with MD simu- 
lations. Our results showed that NPs formed layers which 
started from below the liquid/vapor interface. The struc- 
ture in each layer found in sim ulation s is similar to that 
observed in many experiments j 6 * 10 ^ i.e., a close-packed 
hexagonal 2D crystal. Interestingly, we found that the 
assembly quality was better when the evaporation rate 
was relatively slower. This trend was understood with 
a simple physical picture based on two time scales. The 
time scale associated with NP accumulation in the top 
layer, which is controlled by the evaporation rate, needs 
to be larger than the diffusion time of NPs along the 
interface in order to form a defect-free 2D crystal. Oth- 
erwise, growth of crystals at different locations leads to 
grains with different orientations, and grain boundaries 
and point defects appear when grains meet. Since it is 
known that evaporation needs to fast enough to initiate 
interfacial NP assembly in the first place, our simula- 



tions indicate that there exits an optimum evaporation 
rate at which the packing quality of NP arrays at the liq- 
uid/vapor interface is the best. If evaporation proceeds 
at a higher rate, a relaxation process, which can have a 
large time scale depending on the particle size and in- 
teraction, is needed to "anneal" the crystal and remove 
the imperfections in the packing. Our simulations also 
showed the blockage effect on evaporation because of the 
presence of the dense NP layer near the interface. The ef- 
fective area of the interface was reduced and evaporation 
slowed down compared with a neat solvent. However, in 
our simulations NPs were still immersed in the solvent 
and only the solvent atoms in the region around the top 
of NPs were actually affected in terms of evaporation and 
therefore only a minor reduction in the evaporation rate 
was observed. 
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Supplementary Material 

In this Supplementary Material, we include Voronoi 
constructions of NPs in the top layer for the other 2 
systems that we have studied. They underwent various 
evaporation schemes, as is detailed in the captions of each 
figure. Figure [ST] is for a system with a square cross sec- 
tion (L x = L y = 304a). Figure S2 is for another system 
that had a similar area, but the sizes of the liquid/vapor 
interface were adjusted to accommodate an integer num- 
ber of NPs in both x and y directions (L x = 271a and 



L y = 352a), respectively. 
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(a) (b) (c) 



FIG. SI: Voronoi constructions of NPs in the top layer for the system with L x — L y — 304cr: (a) after 45% (~ 3.2 million 
atoms) of the solvent had evaporated at a fixed rate 2.17 x 10~ 4 r~ 1 a~ 2 for 1.645 x 10 5 r; (b) after 45% (~ 3.2 million atoms) 
of the solvent had evaporated into a vacuum for t = 1.9 X 10 4 r; (c) the same system as in (b) but relaxed for 3.4 x 10 4 r after 
the evaporation was stopped. 




(a) (b) (c) 



FIG. S2: Voronoi constructions of NPs in the top layer for the system with L x — 271a and L y — 352cr: (a) after 42% (~ 3 
million atoms) of the solvent had evaporated at a fixed rate 1.05 x 10~ 3 r~ 1 a~ 2 for 3 x 10 4 r; (b) after 42% (~ 3 million atoms) 
of the solvent had evaporated into a vacuum for 1.6 x 10 4 r; (c) the same system as in (b) but relaxed for 2.8 x 10 4 r after the 
evaporation was stopped. 



